﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.NumericalIntegration
{
    /// <summary>
    /// This class provides an integration for cubic splines.
    /// </summary>
    [Serializable]
    public class CubicSplineIntegrator : AbstractSplineIntegrator
    {
        /// <summary>
        /// This field holds the spline to integrate.
        /// </summary>
        private IFinityCubicSplineIntegral integral;

        /// <summary>
        /// Initializes a new instance of the <see cref="CubicSplineIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The spline to integrate.</param>
        public CubicSplineIntegrator(IFinityCubicSplineIntegral integral)
        {
            this.integral = integral;
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="CubicSplineIntegrator"/> class.
        /// </summary>
        /// <param name="integral">The spline to integrate.</param>
        public CubicSplineIntegrator(FinityCubicSplineIntegral integral)
        {
            if (integral == (FinityCubicSplineIntegral) null)
            {
                throw new ArgumentNullException("integral");
            }

            if (integral.Function == (CubicSpline) null)
            {
                throw new ArgumentNullException("integral");
            }

            this.integral = integral;
        }

        /// <summary>
        /// Gets or sets the spline to integrate.
        /// </summary>
        /// <value>The spline to integrate.</value>
        public IFinityCubicSplineIntegral Integral
        {
            get { return integral; }
            set { integral = value; }
        }


        /// <summary>
        /// Integrates the spline in a numerical way.
        /// </summary>
        /// <returns>The value of the intergral.</returns>
        public override double Integrate()
        {
            bool swapSign = false;
            double a = this.integral.A;
            double b = this.integral.B;

            if (!this.integral.Function.IsDefinedFor(a))
            {
                throw new Exception("The spline is not defined for the integral lower bound a.");
            }

            if (!this.integral.Function.IsDefinedFor(b))
            {
                throw new Exception("The spline is not defined for the integral upper bound b.");
            }

            if (b < a)
            {
                swapSign = true;
                ExtendedMath.Swap(ref a, ref b);
            }

            double result = 0;

            foreach (CubicSplineComponent component in this.integral.Function.Components)
            {
                if ((a <= component.Start.X) && (b >= component.End.X))
                {
                    result += InternalIntegrate(component);
                }
                else
                {
                    if ((component.Start.X < a) && (b < component.End.X))
                    {
                        result += InternalPartIntegrateAToB(component, a, b);
                    }
                    else if ((component.Start.X < a) && (a < component.End.X))
                    {
                        result += InternalPartIntegrateAToEnd(component, a);
                    }
                    else if ((component.Start.X < b) && (b < component.End.X))
                    {
                        result += InternalPartIntegrateStartToB(component, b);
                    }
                }
            }

            if (swapSign)
            {
                return - result;
            }

            return result;
        }

        private static double InternalIntegrate(CubicSplineComponent component)
        {
            SimplePolynomial polynomial = new SimplePolynomial(4);

            polynomial.SetValueAtRank(0, 0);
            polynomial.SetValueAtRank(1, component.Polynomial.ValueAtRank(0));
            polynomial.SetValueAtRank(2, component.Polynomial.ValueAtRank(1) / 2);
            polynomial.SetValueAtRank(3, component.Polynomial.ValueAtRank(2) / 3);
            polynomial.SetValueAtRank(4, component.Polynomial.ValueAtRank(3) / 4);

            return polynomial.SolveAt(Math.Max(component.Start.X, component.End.X)) -
                   polynomial.SolveAt(Math.Min(component.Start.X, component.End.X));
        }

        private static double InternalPartIntegrateAToEnd(CubicSplineComponent component, double a)
        {
            SimplePolynomial polynomial = new SimplePolynomial(4);

            polynomial.SetValueAtRank(0, 0);
            polynomial.SetValueAtRank(1, component.Polynomial.ValueAtRank(0));
            polynomial.SetValueAtRank(2, component.Polynomial.ValueAtRank(1) / 2);
            polynomial.SetValueAtRank(3, component.Polynomial.ValueAtRank(2) / 3);
            polynomial.SetValueAtRank(4, component.Polynomial.ValueAtRank(3) / 4);

            return polynomial.SolveAt(Math.Max(a, component.End.X)) -
                   polynomial.SolveAt(Math.Min(a, component.End.X));
        }

        private static double InternalPartIntegrateStartToB(CubicSplineComponent component, double b)
        {
            SimplePolynomial polynomial = new SimplePolynomial(4);

            polynomial.SetValueAtRank(0, 0);
            polynomial.SetValueAtRank(1, component.Polynomial.ValueAtRank(0));
            polynomial.SetValueAtRank(2, component.Polynomial.ValueAtRank(1) / 2);
            polynomial.SetValueAtRank(3, component.Polynomial.ValueAtRank(2) / 3);
            polynomial.SetValueAtRank(4, component.Polynomial.ValueAtRank(3) / 4);

            return polynomial.SolveAt(Math.Max(component.Start.X, b)) -
                   polynomial.SolveAt(Math.Min(component.Start.X, b));
        }

        private static double InternalPartIntegrateAToB(CubicSplineComponent component, double a, double b)
        {
            SimplePolynomial polynomial = new SimplePolynomial(4);

            polynomial.SetValueAtRank(0, 0);
            polynomial.SetValueAtRank(1, component.Polynomial.ValueAtRank(0));
            polynomial.SetValueAtRank(2, component.Polynomial.ValueAtRank(1) / 2);
            polynomial.SetValueAtRank(3, component.Polynomial.ValueAtRank(2) / 3);
            polynomial.SetValueAtRank(4, component.Polynomial.ValueAtRank(3) / 4);

            return polynomial.SolveAt(Math.Max(a, b)) - polynomial.SolveAt(Math.Min(a, b));
        }
    }
}